##############################################################
#### 0. Load Packages 
##############################################################

packages <- c("ggplot2", "sjlabelled", "labelled", "tidyverse", "cregg", "haven", "dplyr",
              "stargazer", "tidyverse","estimatr","clubSandwich","lmtest", "dotwhisker",
              "scales", "ggthemes", "extrafont","gridExtra")

package_installed <-
  sapply(packages, function(pack)
    pack %in% rownames(installed.packages()))

if (any(!package_installed)) {
  sapply(packages[!package_installed], install.packages)
}

sapply(packages, require, character.only = TRUE)

rm(packages,package_installed)

##############################################################
#### 1.  Load Data 
##############################################################

Sys.setlocale("LC_CTYPE", "en_US.UTF-8")

gep <- load("gep_new.rda")
gep <- labelled::to_factor(gep_new, labelled_only = TRUE, strict = TRUE, unclass = TRUE)

##############################################################
#### 2-1. Conjoint Cleaning Control Group  
##############################################################

process_data <- function(gep, taskNum, choiceNum, controlPolicy) {
  data_prefix <- paste0("c", taskNum)
  control_col1 <- paste0(controlPolicy, "1")
  data <- gep %>%
    dplyr::select(ResponseId, starts_with(data_prefix), starts_with(controlPolicy)) %>%
    mutate(
      taskNum = taskNum,
      choiceNum = choiceNum,
      selected = ifelse(get(control_col1) == paste0("탄소세 정책안 ", LETTERS[choiceNum]), 1, 0),
      redistribution_tax = NA,
      cost = NA,
      energy = NA,
      environment_tax = NA
    )
  
  for (i in 1:4) {
    attrib_name <- paste0(data_prefix, "_attrib", i, "_name")
    attrib_policy <- paste0(data_prefix, "_attrib", i, "_policy", choiceNum)
    
    data <- data %>%
      mutate(
        redistribution_tax = ifelse(get(attrib_name) == "소득 재분배적 목적의 세수 사용", get(attrib_policy), redistribution_tax),
        cost = ifelse(get(attrib_name) == "월 평균 가계 부담<br />(가계 별 월 탄소 1톤 배출 기준)", get(attrib_policy), cost),
        energy = ifelse(get(attrib_name) == "집중 과세 대상 에너지", get(attrib_policy), energy),
        environment_tax = ifelse(get(attrib_name) == "환경적 목적의 세수 사용", get(attrib_policy), environment_tax)
      )
  }
  
  return(data[c("ResponseId", "taskNum", "choiceNum", "selected", "redistribution_tax", "cost", "energy", "environment_tax")])
}

results <- list()
for (task in 1:5) {
  for (choice in 1:2) {
    policyType <- paste0("A", task, "_CQ.Control")
    results[[paste0("dat", task, ".", choice)]] <- process_data(gep, task, choice, policyType)
  }
}
dat_conjoint <- do.call(rbind, results)

dat_conjoint$redistribution_tax <- gsub("<br />", " ", as.character(dat_conjoint$redistribution_tax))
dat_conjoint$cost <- gsub("<br />", " ", as.character(dat_conjoint$cost))
dat_conjoint$energy <- gsub("<br />", " ", as.character(dat_conjoint$energy))
dat_conjoint$environment_tax <- gsub("<br />", " ", as.character(dat_conjoint$environment_tax))

dat_conjoint$redistribution_tax <- as.factor(dat_conjoint$redistribution_tax)
dat_conjoint$cost <- as.factor(dat_conjoint$cost)
dat_conjoint$energy <- as.factor(dat_conjoint$energy)
dat_conjoint$environment_tax <- as.factor(dat_conjoint$environment_tax)

dat_conjoint$cost <- relevel(dat_conjoint$cost, ref = "40,000원")
dat_conjoint$environment_tax <- relevel(dat_conjoint$environment_tax, ref = "환경부 보편 예산에 편입")

dat_conjoint <- subset(dat_conjoint, is.na(selected) != TRUE)


##############################################################
#### 2-2. Conjoint Cleaning Treatment Group
##############################################################

process_data_treat <- function(gep, taskNum, choiceNum, treatmentPolicy) {
  data_prefix <- paste0("c", taskNum)
  treatment_col1 <- paste0(treatmentPolicy, "1")
  treatment_col2 <- paste0(treatmentPolicy, "2_", choiceNum)
  
  data <- gep %>%
    dplyr::select(ResponseId, starts_with(data_prefix), starts_with(treatmentPolicy)) %>%
    mutate(
      taskNum = taskNum,
      choiceNum = choiceNum,
      selected = ifelse(get(treatment_col1) == paste0("탄소세 정책안 ", LETTERS[choiceNum]), 1, 0),
      redistribution_tax = NA,
      cost = NA,
      energy = NA,
      environment_tax = NA
    )
  
  for (i in 1:4) {
    attrib_name <- paste0(data_prefix, "_attrib", i, "_name")
    attrib_policy <- paste0(data_prefix, "_attrib", i, "_policy", choiceNum)
    
    data <- data %>%
      mutate(
        redistribution_tax = ifelse(get(attrib_name) == "소득 재분배적 목적의 세수 사용", get(attrib_policy), redistribution_tax),
        cost = ifelse(get(attrib_name) == "월 평균 가계 부담<br />(가계 별 월 탄소 1톤 배출 기준)", get(attrib_policy), cost),
        energy = ifelse(get(attrib_name) == "집중 과세 대상 에너지", get(attrib_policy), energy),
        environment_tax = ifelse(get(attrib_name) == "환경적 목적의 세수 사용", get(attrib_policy), environment_tax)
      )
  }
  
  return(data[c("ResponseId", "taskNum", "choiceNum", "selected", "redistribution_tax", "cost", "energy", "environment_tax")])
}

treatment_results <- list()
for (task in 1:5) {
  for (choice in 1:2) {
    treatmentType <- paste0("A", task, "_CQ.Treat")
    treatment_results[[paste0("t_dat", task, ".", choice)]] <- process_data_treat(gep, task, choice, treatmentType)
  }
}
t_dat_conjoint <- do.call(rbind, treatment_results)
t_dat_conjoint$redistribution_tax <- gsub("<br />", " ", as.character(t_dat_conjoint$redistribution_tax))
t_dat_conjoint$cost <- gsub("<br />", " ", as.character(t_dat_conjoint$cost))
t_dat_conjoint$energy <- gsub("<br />", " ", as.character(t_dat_conjoint$energy))
t_dat_conjoint$environment_tax <- gsub("<br />", " ", as.character(t_dat_conjoint$environment_tax))

t_dat_conjoint$redistribution_tax <- as.factor(t_dat_conjoint$redistribution_tax)
t_dat_conjoint$cost <- as.factor(t_dat_conjoint$cost)
t_dat_conjoint$energy <- as.factor(t_dat_conjoint$energy)
t_dat_conjoint$environment_tax <- as.factor(t_dat_conjoint$environment_tax)

t_dat_conjoint$cost <- relevel(t_dat_conjoint$cost, ref = "40,000원")
t_dat_conjoint$environment_tax <- relevel(t_dat_conjoint$environment_tax, ref = "환경부 보편 예산에 편입")
t_dat_conjoint <- subset(t_dat_conjoint, is.na(selected) != TRUE)


##############################################################
#### 3. AMCE Estimate                     
##############################################################
#3.1. AMCE Control Group
#3.2. AMCE Treatment Group

### 3.1. AMCE Control Group ##################################

amce.control <- cj(
  formula = selected ~ cost + energy + redistribution_tax + environment_tax,
  data = dat_conjoint, 
  id = ~ ResponseId)

amce.control$level<-c('40,000KRW', '120,000KRW', '80,000KRW', 'Household fuel', 'Industrial fuel', 'Transportation fuel', 'Transfer to General Government Budget', 'Tax Rebate', 'Reducing Corporate tax', 'Reducing Income tax','Subsidy for the Energy-poor', 'Transfer to General Environment Ministry Budget', 'Support Regions Vulnerable to Climate Change', 'Green Job Creation', 'Investment in Invention Renewable Energy')

#AMCE Control Group Binary
plot1<-amce.control %>% 
  ggplot(aes(x=level, y=estimate, color=feature),) + 
  scale_x_discrete(limits=c('40,000KRW', '80,000KRW', '120,000KRW', 'Household fuel', 'Transportation fuel','Industrial fuel', 'Transfer to General Government Budget', 'Tax Rebate', 'Subsidy for the Energy-poor', 'Reducing Income tax','Reducing Corporate tax', 'Transfer to General Environment Ministry Budget', 'Support Regions Vulnerable to Climate Change', 'Investment in Invention Renewable Energy', 'Green Job Creation'))+
  geom_point(position=position_dodge(width = .2)) + 
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  scale_color_manual(name = "Feature of Carbon Tax", values = hue_pal()(4), labels = c("Household Burden(monthly)", "Tax Base Energy Energy", "Redistributive Earmarking","Environmental Earmarking")) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") + coord_flip() + theme(axis.title.x = element_blank()) + theme(axis.title.y = element_blank())+
  ggtitle('Control Group Binary')

### 3.2. AMCE Treatment Group ###############################

t_amce <- cj(
  formula = selected ~ cost + energy + redistribution_tax + environment_tax,
  data = t_dat_conjoint, 
  id = ~ ResponseId)

t_amce$level<-c('40,000KRW', '120,000KRW', '80,000KRW', 
                'Household fuel', 'Industrial fuel', 'Transportation fuel',
                'Transfer to General Government Budget', 'Tax Rebate', 'Reducing Corporate tax', 'Reducing Income tax','Subsidy for the Energy-poor', 
                'Transfer to General Environment Ministry Budget', 'Support Regions Vulnerable to Climate Change', 'Green Job Creation', 'Investment in Invention Renewable Energy')

t_plot1<-t_amce %>% 
  ggplot(aes(x=level, y=estimate, color=feature)) + 
  scale_x_discrete(limits=c('40,000KRW', '80,000KRW', '120,000KRW', 'Household fuel', 'Transportation fuel','Industrial fuel', 'Transfer to General Government Budget', 'Tax Rebate', 'Subsidy for the Energy-poor', 'Reducing Income tax','Reducing Corporate tax', 'Transfer to General Environment Ministry Budget', 'Support Regions Vulnerable to Climate Change', 'Investment in Invention Renewable Energy', 'Green Job Creation'))+
  geom_point(position=position_dodge(width = .2)) + 
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  scale_color_manual(name = "Feature of Carbon Tax", values = hue_pal()(4), labels = c("Household Burden(Monthly)", "Tax Base Energy", "Redistributive Earmarking","Environmental Earmarking")) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") + coord_flip() + theme(axis.title.x = element_blank()) + theme(axis.title.y = element_blank())+
  ggtitle('Treatment Group Binary')



##############################################################
### 4. Plots
##############################################################

amce.control$treatment <- "No Information Treatment"
t_amce$treatment <- "Information Treatment"

amce.all <- rbind(amce.control, t_amce)

pdf("CabonConjointMain.pdf", width=9, height=6)

plot_amce_all<-amce.all %>% 
  ggplot(aes(x=level, y=estimate, color=feature, shape=treatment)) +
  scale_x_discrete(limits=c('40,000KRW', '80,000KRW', '120,000KRW', 'Household fuel', 'Transportation fuel','Industrial fuel', 'Transfer to General Government Budget', 'Tax Rebate', 'Reducing Corporate tax', 'Reducing Income tax', 'Subsidy for the Energy-poor', 'Transfer to General Environment Ministry Budget', 'Support Regions Vulnerable to Climate Change', 'Green Job Creation', 'Investment in Invention Renewable Energy'),
                   labels=c('40,000KRW', '80,000KRW', '120,000KRW', 'Household', 'Transportation','Industry', 'General Government Budget', 'Universial Dividends', 'Corporate Tax Reduction','Income tax Reduction', 'Subsidy for the Energy Poor Population', 'General Environment Ministry Budget', 'Support Regions Vulnerable to Climate Change', 'Green Job Creation', 'Investment in Renewable Energy Development'))+
  geom_point(position=position_dodge(width = .2)) +
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(width = .2)) +
  geom_hline(yintercept=0, linetype="dashed", color="grey") +
  coord_flip() + xlab("Attribute levels") + ylab("Estimate") +
  scale_color_manual(name = "Carbon Tax Attributes", values = hue_pal()(4), 
                     labels = c("Cost for Average Household(monthly)", "Main Tax Target", "Revenue Recycling for Redistribution","Revenue Recycling for Environmental Earmarking")) +
  labs(x = "Attribute levels", y = "Estimate",
       color='Carbon Tax Attributes', shape="Information")+
  guides(shape = guide_legend(reverse = TRUE))+
  guides(color = guide_legend(reverse = TRUE))

plot_amce_all
dev.off()
